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Abstract 



The concepts and methods of Systems Biology are being extended to neuropharmacol- 
ogy, to test and design drugs against neurological and psychiatric disorders. Computa- 
tional modeling by integrating compartmental neural modeling technique and detailed 
kinetic description of pharmacological modulation of transmitter - receptor interaction 
is offered as a method to test the electrophysiological and behavioral effects of putative 
drugs. Even more, an inverse method is suggested as a method for controlling a neural 
system to realize a prescribed temporal pattern. In particular, as an application of the 
proposed new methodology a computational platform is offered to analyze the gener- 
ation and pharmacological modulation of theta rhythm related to anxiety is analyzed 
here in more detail. 
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1 Introduction 



Systems Biology is an emergent movement to combine system level description with 
microscopic details. It might be interpreted as the renaissance of cybernetics Q and 
of system theory |2|, materialized in the works of Robert Rosen |3|. (For an excellent 
review on applying the system theoretical tradition to the new Systems Biology see 


To have a system-level understanding of biological systems 16) we should get 
information from five key features: 

Architecture. The structure (i.e. units and relations among these units) of the system 
from network of gene interactions via cellular networks to the modular architec- 
ture of the brain are the basis of any system level investigations. 

Dynamics. Spatio-temporal patterns (i.e. concentrations of biochemical components, 
cellular activity, global dynamical activities such as measured by electroencephalo- 
gram, EEG) characterize a dynamical system. To describe these patterns dynami- 
cal systems theory offers a conceptual and mathematical framework. Bifurcation 
analysis and sensitivity analysis reveal the qualitative and quantitative changes 
in the behavior of the system. 

Function. This is the role that units (from proteins via genes, cells and cellular net- 
works) play to the functioning of a system (e.g. our body and mind). 

Control. There are internal control mechanisms which maintain the function of the 
system, while external control (such as chemical, electrical or mechanical per- 
turbation) of an impaired system may help to recover its function. 

Design. There are strategies to modify the system architecture and dynamics to get a 
desired behavior at functional level. A desired function may be related to some 
"optimal temporal pattern". 

While Systems Biology is now generally understood in a somewhat restricted way 
for proteins and genes, its conceptual and mathematical framework could be extended 
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to neuroscience, as well. Trivially, there is a direct interaction between molecular and 
mental levels: chemical drugs influence mood and state of consciousness. However, 
"almost all computational models of the mind and brain ignore details about neuro- 
transmitters, hormones, and other molecules." |7|. 

In this paper we show how to realize the program of Systems Biology in the context 
of a new, dynamic neuropharmacology, which was outlined recently elsewhere |8|. 
Also, we offer a methodology to integrate conventional neural models with detailed 
description of neurochemical synaptic transmission in order to develop a new strategy 
for drug discovery. The procedure is illustrated on the problem of finding selective 
anxiolytics. 

In particular, our proposed working hypothesis is that 

• for given putative anxiolytic drugs we can test their effects on the EEG curve by, 

• starting from pharmacodynamic data (i.e. from dose-response plots), which of 
course are different for different recombinant receptors, 

• setting the kinetic scheme and a set of rate constants, 

• simulating the drug-modulated GABA - receptor interaction, 

• calculating the generated postsynaptic current, 

• integrating the results into the network model of the EEG generating mechanism, 

• simulating the emergent global electrical activity, 

• evaluating the result and to decide whether the drug tested has a desirable effect. 

In section |2 we first review the literature and existing computer models of anxiety 
and anxiolytics. Second a review of GABA-GABAa (ligand- receptor) interaction 
is given with a particular attention on its connection to anxiety research. In section [5] 
we elaborate on the proposed working hypothesis and give an example by integrating 
a GABAAreceptor model by Bai et al. (1999)[9 1 into a gamma related theta generating 
interneuron network model by Orban et al. (2001) 1101 . 
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2 Computer models of anxiety 



2.1 Biologically motivated modeling of neural circuits 
2.1.1 Architecture: the septo-hippocampal system 

It was demonstrated (see e.g. the seminal book of Gray and McNaughton 1 1 1 1) that the 
septo-hippocampal system is strongly involved in anxiety and related disorders. 

In a joint pharmacological and computational work 1 12 13 1 effects of the injection 
of the positive and negative GABAa allosteric modulators diazepam and FG-7142, re- 
spectively, were studied. To investigate the dynamical and functional effects of dif- 
ferent pharmacological agents by computational tools a skeleton model of the septo- 
hippocampal system was established. The skeleton network model (Fig.Q of the hip- 
pocampal CA1 region and the septal GABAergic cells consisted of five cell popula- 
tions. The hippocampal CA1 pyramidal cell model was a multicompartmental model 
modified from 1141 . In the hippocampal CA1 region basket neurons and two types of 
horizontal neurons were also taken into account. 

[Figure 1 about here.] 

Connections within and among cell populations were created faithfully following 
the hippocampal structure. The main excitatory input to horizontal neurons is provided 
by the pyramidal cells via AMPA (alpha-amino-3-hydroxy-5-methyl-4-isoxazolepropionic 
acid) mediated synapses 1151 . Synapses of the septally projecting horizontal cells 1161 
and synapses of the other horizontal cell population, the O-LM cell population inner- 
vating distal apical dendrites of pyramidal cells 11171 are of the GABAa type. O-LM 
neurons (on of the two horizontal interneurons) also innervate parvalbumin contain- 
ing basket neurons HI 81 . Basket neurons innervate pyramidal cells at their somatic 
region and other basket neurons 1 19 1 as well. Septal GABAergic cells innervate other 
septal GABAergic cells and hippocampal interneurons [20 21 1 (Fig. 0. For a full 
description of this model see the online supplementary materials to the paper 111 21 at: 
http://geza.kzoo.edu/theta/theta. html. 

The above described model captures several elements of the complex structure of 
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the hippocampal CA1 and can be used to account for very precise interactions within 
this region. However, when the focus of interest is rather on general phenomena tak- 
ing place during rhythm generation modelers might settle for a simpler architecture. 
In 1 10 1 authors describe gamma related theta oscillation generation in the CA3 region 
of the hippocampus. The architecture of the model is exceedingly simplified: only an 
interneuron network is simulated in detail. This simplification, however, allowed the 
authors to introduce an extrahippocampal input and study its effect on rhythm gener- 
ation. As a result, the model is able to account for basic phenomena necessary for 
the generation of gamma related theta oscillation. As an extension of this model, au- 
thors show [ 22 1 that activity of the interneuron network indeed accounts for rhythm 
generation in pyramidal neurons. 

2.1.2 Dynamics: Generation of theta rhythms 

Theta frequency oscillation of the septo-hippocampal system has been considered as a 
prominent activity associated with cognitive function and affective processes. To inves- 
tigate the generation and control of theta oscillation in the hippocampal CA1 region the 
previously described detailed, realistic model 1121 was used. Firing of neurons of the 
four simulated populations were modulated in time. There are time intervals in which 
firing was significantly reduced were alternated by intervals where enhanced firing was 
observed. This synchronized state of neural firing was further confirmed by the field 
potential, which exhibited a prominent «5 Hz oscillation (see Fig. 6 in 1 12 1). 

Simulation results showed that key components in the regulation of the popula- 
tion theta frequency are the membrane potential oscillation frequency of pyramidal 
cells, strength of the pyramidal cell - O-LM cell innervation and the strength of recur- 
rent basket cell connections. Amplitude and frequency of pyramidal cell membrane 
potential oscillation is determined by their average depolarization, passive membrane 
parameters and parameters of the active currents. Average depolarization in our model 
results from septal cholinerg innervation. An important factor is the presence and max- 
imal conductance of the hyperpolarization activated current. If Ih is present it shortens 
response times of pyramidal cells to hyperpolarizing current pulses and more impor- 
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tantly decreases its variance: Ij, acts as a frequency stabilizer. Synaptic strengths in our 
description are set by convergence numbers and maximal synaptic conductances. 

An explanation of intrahippocampal theta oscillation generation - based on this 
model - includes (i), signal propagation in the pyramidal cell — > O-LM cell — > basket 
cell — > pyramidal cell feed-back loop, (ii), synchronization of neural activity via the 
recurrent, inhibitory GABAa connections within the basket cell network and (iii), syn- 
chronization of pyramidal cell firing due to rebound action potential generation. It is 
true that the propagation of a single signal throughout this trisynaptic loop would not 
require the amount of time characteristic to the theta oscillation (wO.2-0.25 sec), thus 
in the present case the population oscillation is created not by the propagation of single 
signals but rather the propagation of a "synchronized state" in the network. 

The observed periodic population activity is brought about by alternating synchro- 
nization and desynchronization of cell activities due to the interplay of the above men- 
tioned synchronizing forces and some desynchronizing forces (such as heterogeneity 
of cell parameters and diversity of synaptic connections), as observed in 11 01 1221 . 

2.1.3 Function: Mood Regulation 

The hippocampal formation is known to be involved in cognitive processes (navigation 
and memory formation - for reviews see e.g. l23i 24 1) and mood regulation | lTI 1251 . 
When an impairment occurs in the architecture or dynamics of the septo-hippocampal 
system, this change is reflected in its function e.g. as an impairment of its memory 
storing or mood regulatory capability. One such mood disorder is anxiety: 

'Anxiety is a complex combination of the feeling of fear, apprehension and worry 
often accompanied by physical sensations such as palpitations, chest pain and/or short- 
ness of breath. It may exist as a primary brain disorder or may be associated with other 
medical problems including other psychiatric disorders. . . 

A chronically recurring case of anxiety that has a serious affect on your life may be 
clinically diagnosed as an anxiety disorder. The most common are Generalized anx- 
iety disorder, Panic disorder, Social anxiety disorder, phobias, Obsessive-compulsive 
disorder, and post-traumatic stress disorder" ll26l 
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It is well documented that anxiolytics and hypnotics reduce the amplitude of septo- 
hippocampal oscillatory theta activity, which contributes to their therapeutic effect but 
causes unwanted side effects, e.g. cognitive impairment as well 1 27 1 1281 . 

Historically used mood regulators act on the barbiturate or benzodiazepine (BZ) 
sites of GABA receptors. These drugs have both anxiolytic and hypnotic activity be- 
cause they modulate all types of GABAa receptors. 

Conventional anti-anxiety drugs have sedative-hypnotic side effects, and have neg- 
ative effects on cognitive functions, such as memory and generally they reduce the am- 
plitude of theta rhythm. While it seems possible to dissociate sedative-hypnotic effects 
and anxiety, it is not clear how to decrease anxiety without impairing memory function. 

Models of anxioselective actions: search for data 

Recently it became clear that a subunits of GABAa receptors exhibit a remarkable 
functional specificity. Genetic manipulations helped to show that oci subunits are re- 
sponsible for mediating sedative effects, while 0C2 subunits mediates anxiolytic effects 
1291 . Preliminary experimental data and modelling studies for the effects of the pref- 
erential GABAa cci and 0.2 positive allosteric modulator, Zolpidem and L838,417, re- 
spectively, for the septo-hippocampal theta activity have been reported 1 30 1 . 

In that study we examined the effects of the cti and 0C2 subtype-selective BZ site 
ligand Zolpidem and L838,417 on the septo-hippocampal system. In electrophysiolog- 
ical experiments extracellular single unit recordings were performed from medial sep- 
tum/diagonal band of Broca with simultaneous hippocampal (CA1) electroencephalo- 
gram (EEG) recordings from anesthetized rats. Both of the drugs eliminated the hip- 
pocampal theta oscillation, and turned the firing pattern of medial septal cells from 
periodic to aperiodic, but only the Zolpidem reduced the firing rate of the these neu- 
rons. In parallel to these experimental observations, a computational model has been 
constructed to clearly understand the effect of these drugs on the medial septal pace- 
maker cells. We showed that the aperiodic firing of hippocampo-septal neurons can 
reduce the periodicity of the medial-septal cells, as we have seen in the case of the 
L838,417. The reduction of firing rates in the case of Zolpidem is attributed to the 
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increase of the synaptic conductances and the constant inhibition of these cells. We 
modelled these drug effects by modifying (i) the synaptic maximal conductances of the 
GABA synapses, (ii) the constant excitatory drive of the median septal cells and (iii) 
the hippocampal input. The incorporation of a more detailed synaptic model, similarly 
to the one proposed in Section l3~Tl taking into account differences between oti and 0C2 
subunits is in progress. 

Zolpidem increases by concentration-dependent manner the duration and amplitude 
of the postsynaptic current, most likely by enhancing the affinity of the receptors for 
GABA 1 3 1 1, but these effects were diminished or absent in neurons from cti knock-out 
mice 1 32 1 

There seem to be compounds, which might have comparable binding affinity but 
different efficacies at the various subtypes, thereby preferentially exerting their effects 
at subtypes thought to be associated with anxiety. L838,417 seems to be an example for 
subtype selective compounds 1 33 1, but only few kinetic or even pharmacodynamic data 
could be found in the public domain. In Section[3]drugs whose effect on the GABA ki- 
netics is known will be used to demonstrate the proposed methodology, using a detailed 
kinetic model of GABA receptors within the network model of CA3 interneurons. 

2.2 Phenomenological and kinetic descriptions of synaptic interac- 
tion 

2.2.1 Control: GABA receptor kinetics 
Pharmacological modulation: why and how? 

One perspective of the pharmacotherapy of anxiety might be based on reducing septo- 
hippocampal theta activity and it is related to modulate the interaction between the 
GABA transmitter and GABAa receptors. Since GABAa receptors occur in a variety 
of specific forms, they might be subject of pharmacological control to obtain selective 
functional effects. While by-and-large it is understandable that putative drugs induce 
opening and closing of synaptic channels (here the chloride channels), the underlying 
finely tuned kinetic mechanism of the interaction of transmitters, modulatory agents 
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and receptors are not known. Our longer term goal is to find out how pharmacologi- 
cal effects depend on the concentration of transmitter and of modulators. In technical 
terms, while there are detailed models of the presynaptic transmitter release here we 
restrict ourselves to build a detailed kinetic model for the generation of the postsynap- 
tic response. 

Pharmacological elements 

Receptor Structure: GABAa receptors are pentameric structures consisting of multi- 
ple subunits. At this moment 1341 nineteen subunits have been cloned from mammalian 
brain. According to their sequence similarities, subunits have been grouped into seven 
families: a, (3, y, 6, e, n and 9. Only a few dozen among the many combinatorial pos- 
sibilities exist. The most frequent subtypes are two ot, two (3 and one y subunits. The 
structural variations imply functional consequences |34|, among others for the kinetic 
properties. 

Pharmacodynamics: Effect of drugs can be characterized by dose-response curves. 
Some drugs have the ability to open a certain ion channel per se. In this case drugs are 
characterized by the maximal evoked response or the potency, the drug concentration, 
which elicit 50 % of the maximal response. 

There are other cases, however, when the drug can not open the ion channel directly 
but requires an additional ligand to evoke an action. In the followings we study drugs 
that modulate interaction between the GABA transmitter and GABAa receptors. 

Both barbiturates and BZs shift the GABA concentration-response curve to the left, 
but barbiturates also increase the maximum response. They act on different states, con- 
sequently they have different kinetic effects: average open time of the channel, but not 
the channel opening frequency is increased significantly by barbiturates. As opposed 
to BZs, barbiturate receptors do not contain y subunits (see later). One more differ- 
ence is that at high concentration GABA receptor channels can directly be opened by 
barbiturates. For a summary see 1351 . Anxiolytic activity was not a particular disad- 
vantage when these drugs were used as hypnotics, hypnosis was a definite disadvantage 
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when they were used as anxiolytics. Recent discoveries made possible the separation 
between hypnotic and anxyolitic activity and selective hypnotic agents (e.g. Zolpidem) 
are already on the market. Selective anxiolytics are on the preclinical and/or in clinical 
trial stage. 

The conventional tool of computational neuroscience: a phenomenological model 

One way to describe synaptic transmission is to use a gating variable similar to the well 
known Hodgkin - Huxley formalism: 



Isyn = g S ynS(V-E syn ) (la) 

ds 

_ = cxF(V pre )(1-s)-(3s (lb) 
r (V pre ) = J — 5 x (lc) 

V nre Cs' 



1 + exp 



' pre ^syn 



K 



with I syn being the synaptic current, g syn the maximal synaptic conductance, s the 
gating variable of the synaptic channel, E syn the synaptic reversal potential, F(-) is an 
activation function, a and |3 rate functions describing opening and closing of the gate 
of the synaptic channel, syn is a threshold. 

Figure [2] illustrates the general form of effects of GABAa receptor modulators 
described by the phenomenological model, where the effect of drugs was taken into 
account solely by modifying g syn , the maximal synaptic conductance. 

[Figure 2 about here.] 

From pharmacodynamics to detailed kinetic schemes 

Kinetic models of transmitter-receptor interaction (without and in the presence of drugs) 
help to calculate the induced synaptic current, and gives a possibility to get quantita- 
tive insight how to modulate it. A specific step in the transmitter-receptor interaction, 
namely desensitization, plays an important role in shaping the GABAergic currents. 
Since our general goal is to incorporate appropriate kinetic models to the network 
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model to compute the drug effects on global brain rhythm, we should have appropriate 
kinetic model. 

As opposed to the conventional phenomenological description, a more effective, 
but certainly most expensive, modelling tool to evaluate the pharmacological effects 
of the different modulators, or even to give help for offering new putative molecules 
for drug discovery, is the inclusion of more detailed kinetic studies of GABA receptor 
modulation. 

Jones and Westbrook 1361 established a model for describing the rapid desensiti- 
zation of the GAB A a receptors. More specific kinetic models should be studied to 
describe the effects of the different (full and partial) agonists and antagonists taking 
into account that connections between singly and doubly bound open and desensitized 
states may influence the receptor kinetics. Baker et al. |37| explained the functional 
difference between the effects of propofol (which has hypnotic effect) and of midazo- 
lam (a sedative - amnestic drug) based on a detailed kinetic model of a single cell with 
autapses and two interconnected cells. 

The main difference is that propofol modifies the desensitization processes, more 
dramatically the slow desensitization steps and the modified kinetic parameters. These 
differences imply distinct behavior of the network (synchronization, frequency of os- 
cillation) and therefore also in function. 

3 Integrating receptor kinetics and neural network mod- 
els 

3.1 Setting up the model 

In the following section we connect the two concepts reviewed in the first part of this 
paper. Below is the description of an illustrative model integrating the detailed kinetic 
description of synaptic receptors |9| into the biophysical model of a gamma related 
theta rhythm generating interneuron network 1101 1221 . We use this model to study the 
effect of drugs that have a well identified effect on the chemical reactions taking place 
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in the GABAa receptor. 

Here, the model is introduced only briefly, detailed equations and parameters can 
be found in the Appendix. The modeled interneurons contained the Hodgkin - Huxley 
channels (sodium, potassium and leak), which were identical for all cells (see lA-lal - 
lA-lmt . Also, the synaptic current (I syn (t)) and an input current (1^ (t)) were considered 
(see below). 

A connection between two interneurons was established through a synaptic mecha- 
nism consisting of a phenomenological presynaptic and a chemically realistic described 
postsynaptic part. The presynaptic model describes the transmitter release due to an ac- 
tion potential generated in the presynaptic neuron by a sigmoid transfer function (see 
IA-2> . The value of the released transmitter concentration ([L]) is then used in the post- 
synaptic GABAa receptor model to calculate the concentration of receptor proteins 
being in the conducting (two ligand binded open) state ([L20](t)) (see IA-31 — |A-4g^ . 
The synaptic current was thus given by 

N 

Isy„(t) = Y_ 9syn,i • [L20] t (t) • (V(t) - V syn ) (2) 

1=1 

where 1M is the number of presynaptic neurons, i indexes these presynaptic neurons, 
g syn is the maximal synaptic conductance representing the synaptic channel density, 
V(t) is the postsynaptic membrane potential and V syn = —75 mV is the reversal po- 
tential of the chloride current. For details and parameters of the postsynaptic model see 
the Appendix. 

Using this synaptic mechanism interneurons were connected into a random net- 
work. Probability of forming a connection between any two neurons was from 60 % 
to 100 % in the simulations; autapses were not allowed. Maximal synaptic strength was 
normalized to enable scalability of the model: for each postsynaptic neuron Y.i=] 9syn,i = 9s 
where g syn was varied from 0. 1 mS/cm 2 to 0.2 mS/cm 2 . 

As it was shown in previous computer models 1 10| 1221 when the interneurons are 
driven by periodically modulated input and the phase of inputs to different cells is 
heterogeneous, the random network of interneurons generate gamma related theta fre- 
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quency resonance. Indeed, as shown before by Wang and Buzsaki (1996) a network 
of mutually interconnected fast firing inhibitory neurons can generate synchronized 
gamma frequency oscillation. This synchronization is reflected e.g. in the population 
activity or in the raster diagrams (for these measures see the Appendix). However, 
when cells are driven by a periodic inputs with disperse phase instead of a constant 
current, a slow (theta frequency range) modulation of the population activity emerges 
(for details on necessary and sufficient conditions for theta generation the reader is 
kindly directed to Orban et al. (2001)). Briefly, the input for each cell took the follow- 
ing form: 

Ii(t) =DC-sin(cut + 4K) + DC (3) 

with DC being the amplitude and the offset of this current (0.005 - 0.01 A/m 2 ), w the 
frequency (25 - 50 Hz) and 4h is the phase drawn from a Gaussian distribution with 
mean and 30 degree variance. 

3.2 Direct problem: in silico drug screening with the integrated 
model 

Kinetic modeling of synaptic transmission has a flexibility in the level of detailed de- 
scription from chemical kinetic to simplified representation 1 38 1. The development of 
new pharmacological, electrophysiological and computational techniques make possi- 
ble to investigate the modulatory effects of putative drugs for synaptic currents, and 
consequently for local field potentials and even behavioral states. Putative drugs with 
given kinetic properties can be tested in silico before (instead of?) real chemical and 
biological studies. 

As noted before population oscillations reflected by the EEG can be considered as 
biomarkers for some brain disorders. We used our integrated model to study effects of 
two known drugs, the hypnotic propofol and the sedative midazolam on the population 
theta rhythm. In a previous work 1221 it was shown that the population activity of the 
interneuron network can be used to account for the activity of pyramidal cells, which 
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can be in turn translated into EEG. Thus, to quantify our results the raster plots, a 
population activity measure and its Fast Fourier Transform were used (see Appendix). 

Control condition: As shown previously 1 10 22 1 in control situation this model 
exhibits spike synchronization in the gamma frequency band, which is modulated in 
the theta band (Fig. [3} (for a detailed explanation of measures used below see the 
Appendix). On Fig.|3jl (the raster plot) each row symbolizes the spike train of a given 
cell with a black dot representing the firing of an action potential. Note that firing 
of the cells line up in well-synchronized columns e.g. from approximately 2.2 sec to 
2.35 sec. As a result of heterogeneity in the phase of the driving current and in synaptic 
contacts this high synchrony loosens, column on the raster plot are less well organized 
in the 2.35 sec-2.45 sec interval. This synchronization -desynchronization of cells is 
a result of the modulation of the instantaneous frequency of individual cells (Fig. |5J? 
- for more details the reader is kindly referred to Orban et al (2001)). A qualitative 
measure of the synchrony is given by the population activity (Fig.|5Jr) and the network 
synchrony (Fig.|3p). 

[Figure 3 about here.] 

As argued in \22 \ this synchronization- desynchronization of the population activ- 
ity of interneurons is transposed to principal cells of the hippocampal CA3 modulating 
their firing probability in the theta band. This modulation might be reflected in the 
EEG or CA3 local field potential. Thus, we propose that there is a mapping between 
the FFT of the CA3 local field potential or EEG and the FFT of the simulated popula- 
tion activity (Fig.|3fi), which shows a peak at the theta frequency. Other peaks at higher 
frequencies are related to the fast firing properties and the resulting synchronization of 
the interneurons as well as the frequency of the driving current. 

Drug administration: Given the detailed kinetic description of GABAa receptors 
the modulation of the synaptic current due to its interaction with the drug can be simu- 
lated. Moreover, in the integrated model framework an evaluation of the drug effect on 
the network level can be performed. 

The identification of the effect of a given drug on the kinetic rate parameters of a 
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certain receptor type or subtype is a complicated and tedious undertaking. To present 
the soundness of our model two known drugs were used in the simulations: propofol 
and midazolam (for parameter values see the Appendix). However, a considerable ef- 
fort was made to estimate the change of the rate constants due to the application of 
the L838,417 compound using the chemical method of parameter estimation. Unfor- 
tunately, sufficient kinetic or even pharmacodynamic data could not be found in the 
public domain to obtain realistic estimation for the rate constants. 

Figure |4] and |5] show effects of propofol and midazolam, respectively, on the net- 
work behavior (sub-figures are the same as in Fig.[5J A: raster plot; B: the instantaneous 
frequency versus time; C: the population activity versus time; D: network synchrony; 
E: FFT of the population activity function). The purpose of the present paper is not to 
give a detailed description and analysis of the mode of action of these drugs, thus we re- 
strict ourselves to mention that while propofol prolongs deactivation of the receptor and 
reduces the development of both fast and slow desensitization, midazolam facilitates 
the accumulation of the receptors into the slow desensitized state, which compensates 
the current resulting from slower deactivation |9|. As a result propofol enhances the 
tonic synaptic current to a greater extent than midazolam, whereas propofol and mi- 
dazolam produces similar changes to the time course of single inhibitory postsynaptic 
currents (IPSCs) 1591 . 

[Figure 4 about here.] 
[Figure 5 about here.] 

On the network level the difference that propofol and midazolam exert on the tonic 
component of the synaptic current has a grave consequence. While the peak in the FFT 
at approximately 4 Hz representing the theta modulation of the population activity is 
still present in the case when synapses were "treated" with midazolam (Fig. |5ji) it is 
completely missing for the propofol "treated" case (Fig. 

These findings are in concert with our earlier results |12| where we have shown 
that diazepam a positive allosteric modulator of the GABAa synapses inhibited theta 
oscillation. There synapses were modeled using the phenomenological technique and 
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the increase of the tonic component of the synaptic current was achieved by increasing 
the maximal synaptic conductance parameter, while the time constants were kept con- 
stant to account for similar IPSCs. Methodologically, the approach we followed in the 
integrated model is more established, as generally the maximal synaptic conductance 
parameter is mostly associated with channel density in the synapse, which is not mod- 
ified by the administration of the drug. Moreover, using the detailed kinetic model of 
the GABAa synapse allows us to i, account for drug induced changes in the receptor in 
a chemically more sound way; ii, simulate a much larger variety of situations than by 
using the phenomenological model; and Hi, present several possibilities (several sets of 
rate constants) resulting in the same receptor behavior (a given receptor behavior might 
be brought about by a non-unique set of parameters) to drug researchers who in turn 
can engineer the desired molecule that is easiest to design. 

4 Discussion 

4.1 Design (Inverse Problem): From System Identification to Op- 
timal Temporal Patterns 

We have shown that in a moderately complex conductance-based model of the hip- 
pocampal CA3 region theta rhythm generation can be observed and its general proper- 
ties can be identified 1 10 1 These results qualify the model for consideration as a useful 
tool in the hands of pharmacologists, physiologists and computational neuroscientists 
to complete their repertoire of available tools in the search for efficient and specific 
drugs. 

[Figure 6 about here.] 

Figure|6]is an oversimplified scheme offered for finding a modulator to set optimal 
septo-hippocampal EEG pattern. 

In order to decrease anxiety first a desired EEG pattern should be defined. Anxy- 
olitics should reduce the theta amplitude (but preserve the cognitive performance and 
avoiding sedative and hypnotic side effects). Computational analysis should offer a 
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kinetic scheme and the best rate constants to modulate the fixed network to minimize 
the deviation from the desired "optimal pattern". (Network architecture is supposed to 
be fixed. By neglecting this assumption we should turn from neuropharmacology to 
neurosurgery...) Most likely there are more than one possibilities to reach the goal, and 
model discrimination and parameter estimation techniques may help to narrowing the 
alternatives. 

As it is known from chemical kinetics 1401 14 II sensitivity analysis shows that in a 
kinetic scheme there are "more and less important" components and reactions. It may 
help to answer the question how to modify a given drug to change the reaction rate 
constants in the desired direction - and leaving everything else intact. 

4.2 Further research 

The aim of the present paper is to offer conceptual and mathematical frameworks to 
integrate network and receptor level descriptions for investigating the effects of po- 
tential drugs on the global electrical patterns of a neural center, and on the behavioral 
states (mood, consciousness, etc.). Once we have understood (i) the basic mechanisms 
of rhythm generation, (ii) the elementary steps of the modulatory process, we shall be 
able to give advice to drug designers pointing out which subprocess and how to be 
modulated to reach the given goal. 

Specifically, we briefly reviewed some aspects of GABAa receptor kinetics, and 
the effects of (full and partial) agonists, antagonists and inverse antagonists to septo- 
hippocampal theta rhythms. Our specific goal is to offer a computational platform to 
design anxiolytic drugs with as small as possible side effects. While it is known that 
positive allosteric modulators acting on GABAa <Xi subunits are potential candidates 
for being selective anxyolitics, integrative computational modeling would help to select 
potential drugs with the appropriate kinetic properties. 
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APPENDIX 



Interneuron model 

The interneuron model was taken from |42| and obeys the following current balance 
equation: 



dV 

C m — = -I Na - Ik - t - V + It (A-la) 

lL = g L (V-E L ) (A-lb) 

Ine = gN a mooh.(V- E Na ) (A-lc) 

TTtoo = — — (A-ld) 

^-m Pm 

_ -O.KV + 35) 

m " exp(-0.1 (V + 35))-1 ( ' 

13,11 = exp (-0.1 (V + 28)) (A " lf) 



dh 
dt 

a h = 0.07 exp ( (V J~ 58) J (A-lh) 



4> («h (1 — K) — fihh) (A-lg) 
20 

exp (—0.1 (V + 28)) 

lK = g K n 4 (V-E K ) (A-lj) 
dn 
dt 



cMcUl -n)-P„n) (A-lk) 



-O.KV + 34) 
an " exp (-0.1 (V + 34))-1 ( ' 

(3 n ^0.125exp ^ (V 8 + 44) ") (A-lm) 



with parameters: g L = 1 S/m 2 , g Na = 350 Sim 2 , g K = 90S/m 2 , E L = — 65.3mV, 
E Na = 55 mV, E K = -90 mV cb =5, 



20 



Synapse model 

The presynaptic model describes transmitter release due to presynaptic action poten- 
tials. The following sigmoid function was used: 

where [L] (t) is the released GABA, V pre (t) is the presynaptic membrane potential. 

Our starting point for the postsynaptic model is the work of Bai et al. 1 9 1 originally 
developed by Celentano and Wong [43 1 to describe the emergence of open channels as 
the effect of GABA binding to receptors. However, we explicitly show the presence of 
an essential participant of the model, the ligand L differently from in Scheme 1 of the 
cited paper and our previous work |40|. Obviously, in the model the ligand L denotes 
GABA, whereas I_2 O denotes the open channels, the notation expressing the fact that 
two ligand molecules are needed to change the receptor into an open channel: 



L 2 D 



fast 

d f % r f 

2ko n k on d s 

L + C ^ LiC L + L,C ^ L 2 C ^ L 2 D slow (A-3) 
Ks 2koff r s 

a UP 

L 2 
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The induced differential equations are the followings: 



[C] ' = -2ko n [L] [C] + KdU C] (A-4a) 

[Li C] ' = 2k on [L] [C] + 2k off [L 2 C] - (k off + k on [L] ) [Li C] (A-4b) 

[L 2 C]' = ko [L][Li C] +r f [L 2 D fast ] +r s [L 2 D slow ] + a[L 2 0] (A-4c) 

- (2k off + d f + d s + |3 ) [L 2 C] (A-4d) 

[L 2 D fast ] ' = d f [L 2 C] - r f [L 2 D fast ] (A-4e) 

[L 2 D slow ] ' = d s [L 2 C] - r s [L 2 D slow ] (A-4f) 

[L 2 0]' = (3[L 2 C]-a[L 2 0]. (A-4g) 

Furthermore, we have the following conservation equation for the total quantity of 
channels in any form: 

[C] + [Li C] + [L 2 C] + [L 2 D fast ] + [L 2 D slow ] + [L 2 0] = [C](0) = 1 . (A-5) 



Rate constants in the model were set up according to 1 37 1 to account for the control 
situation, the propofol treatment and the midazolam treatment: 
Rate Constants (ms -1 ) 



Parameters Control Propofol Midazolam 



kon 


1000/M 


1000/M 


1000/M 


k ff 


0.103 


0.056 


0.056 


d» 


3.0 


1.62 


3.0 


rf 


0.2 


0.12 


0.2 


d s 


0.026 


0.014 


0.026 


r s 


0.0001 


0.0001 


0.0001 


a 


0.4 


0.4 


0.4 


P 


6.0 


6.0 


6.0 



Data analysis 
The raster plot 

A raster plot is used to visualize spike-trains of neurons. To draw a raster plot, first 
spiking times of neurons were identified, then time was discretized assigning a or 1 
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value to each time bin depending on weather the cell emitted a spike in a given bin or 
not: 



t 



spiking 



{t | Vt(t) >0AV t (t) >0} 



F7(l) = 



t 



spiking _ f 



n{t | It < t < (1 +1)t} 



(A-6a) 
(A-6b) 



where t^ plking is the set of time points when cell i emits an action potential, Vj.(t) is the 
membrane potential of cell i, F^(l) is the spike train and t = 5 ms is the binning width. 



The instantaneous frequency 

To characterize the behavior of individual cells a measure of their instantaneous fre- 
quency or firing rate was calculated the difference as in 1 10 1: 



ISIt (t) = min (tf ikmg D t 2 |t 2 > t) - max (if** n U |ti < t) (A-7a) 

fi(t) = — (A-7b) 
1 ' ISIt(t) v ; 



The population activity measure 

A "projection" of the raster plot to the time axis gives the population activity: 

a T ' T '(t) = -^ Y_ (A-8) 
1=1 i=[|] 

where N is the number of neurons, T' = 10 ms is a time window in which action 
potentials are counted and [m] denotes the integer part of m. 
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The network synchrony 

The network synchrony measure was defined as follows: 



N v-1 



K T ' T 'm 



N(N - 1 



^11 



i=2 j = 1 



(A-9) 



\ i=[*] i=[f] 



where the window width was set to T' = 50 ms. 
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